XY Spin Fluid in an External Magnetic Field 
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A method of integral equations is developed to study inhomogeneous fluids with planar spins in an 
external field. As a result, the calculations for these systems appear to be no more difficult than those 
for ordinary homogeneous liquids. The approach proposed is applied to the ferromagnetic XY spin 
fluid in a magnetic field using a soft mean spherical closure and the Born-Green- Yvon equation. 
This provides an accurate reproduction of the complicated phase diagram behavior obtained by 
cumbersome Gibbs ensemble simulation and multiple histogram reweighting techniques. 
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Spin fluids are examples of many body systems show- 
ing a rich variety of phases in the global phase diagram 
[1-4]. Besides gas-liquid (G-L) and para-ferro-magnetic 
(P-F) phase transitions, tricritical, critical end and triple 
point behavior is observed. Under special conditions, 
an unsymmetrical tricritical van Laar point exists ad- 
ditionally [5]. This complexity arises due to a coupling 
between spin and spatial interactions. Similar phase di- 
agrams are found in symmetric binary mixtures [6-10] 
with their deniixing and G-L transitions, spin lattice gas 
models [11,12], mixtures of '^He-^'He with the superfluid 
and demixing states [13-15], and other systems. 

The properties of spin fluids were studied using mean 
field (MF) theories [1-4] , more accurate integral equation 
(IE) approaches [5,16-20], and Monte Carlo (MC) simu- 
lation techniques [4,16,19,21-24]. Different types of mod- 
els, such as the well-recognized discrete ID spin Ising, or 
continuous 2D spin XY and 3D Heisenberg fluids, have 
been considered. Despite this, the question concerning 
the global phase diagram topology of the XY spin fluid 
including the influence of an external magnetic field has 
never been addressed. Moreover, the IE approach has 
been restricted either to simplified ideal Heisenberg fluids 
[16-20], where nonmagnetic attractive interactions are 
absent, or to ideal and nonideal Ising models [5]. 

Surprisingly, up to now there were no attempts on de- 
veloping the IE approach for the XY spin fluid model. 
This model may play a crucial role in the description 
of superfluid transitions in pure ^He and its mixtures in 
bulk or in media such as porous gold [15] or silica aero- 
gel [25]. It is generally believed [25] that the superfluid 
transition in "'He belongs to the classical 3D XY model 
universality class (here 3D relates to the dimensionality 
of spatial coordinates). On the other hand, the fluid of 
particles with embedded XY spins described by classical 
statistical mechanics can be treated as one of the sim- 
plest model of disordered continuum systems exhibiting 
ferromagnetic behavior. 

The presence of spin interactions and external fields 
destroys the homogeneity of the fluid, producing nonuni- 
formity or anisotropy in the one-body density. Within 
the standard IE approach this leads to the necessity of 



performing very complicated joint calculations for one- 
and two-body distribution functions on the basis of the 
coupled set of the inhomogeneous Ornstein-Zernike (lOZ) 
equation, a closure relation, and the first equation of the 
Born-Green- Yvon (BGY) hierarchy [26]. Such calcula- 
tions result in unresolvahle numerical difficulties because 
of the restricted capabilities of modern supercomputers. 
It is worth emphasizing also that existing IE develop- 
ments for Ising [5] and Heisenberg [16-20] systems are not 
applicable to the XY fluid. The reason is that neither it 
can be mapped onto a binary homogeneous mixture (as 
for Ising) nor its anisotropic correlations be expanded in 
spherical harmonics (as for Heisenberg). The specific XY 
spin interactions require a separate IE consideration. 

In this Letter we present a method allowing to over- 
come the difficulties of the lOZ approach in the case of 
XY fluids. Comparison of the obtained IE solutions with 
our simulation results has shown a quantitative reproduc- 
tion of the phase diagrams in a wide region of tempera- 
ture, density, external field and interaction parameters. 

Consider an XY spin fluid model with the Hamiltonian 
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where iV is the total number of particles, is the 3D 
spatial coordinate of the i-th body carrying 2D spin 
of unit length, Vij = [r^ — r^l denotes the interparticle 
separation, and H is the external magnetic field vector 
lying like in the XY-plane. The exchange integral J of 
ferromagnetic interactions and the nonmagnetic attrac- 
tion potential / can be chosen in the form of Yukawa 
functions, 

J(r-) = -exp(-^), I{r)^J{r)/R, (2) 

where e and a denote the interaction intensity and the 
size of the particles, respectively, with R being the ratio 
defining the relative strength of J to /. The repulsion 
(j) between particles can be modeled by a more realistic 
soft-core (shifted Lennard- Jones) potential [4,5], 
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r < ^a, 
r > v^cr, 
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rather than by the hard sphere one. 

A complete thermodynamic and magnetic description 

of system (1) can be performed in terms of orientationally 
dependent one-body £,{tp) and two-body g{r,(pi,(p2) = 
h{r, (fi,(p2) + 1 distribution functions. The angles are 
referred to the external field, so that H • s = _ff cos (p and 
Si- S2 — cos{ipi — 1P2). According to the liquid state the- 
ory [26], the total correlation function h satisfies the lOZ 
equation which in our case reads 
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where p = N/V is the particle number density, V the 
volume and c{r,(fii,(p2) the direct correlation function. 

The lOZ equation (4) must be complemented by a clo- 
sure relation. The most general form of it is 



g = exp( — l3u + h — c + 



(5) 



where u{r, tpi , (^2) = (pir) — I{r) — J{r) cos(<pi — ip2) with 
/3~^ = k^T being the temperature, and B is the bridge 
function. This function cannot be determined exactly 
for any system of interacting particles, but a lot of ap- 
proaches exist allowing to present it approximately [26] . 
One way is to use the soft mean spherical approximation 
(SMSA) [5,27] 

B{r,Lpi,ip2) = ln[l -|-r(r,(/Ji,(^2)] - T(r, (^1, (/?2) , (6) 

where t = h — c — pu\. The long-ranged part u\ can be 
extracted [5] from the total potential u as u\{r, 1^1, (^2 ) = 
-[/(r) +J{r) cos((yJi - (^2)] exp[-/?0(r)]. 

Evaluation of pair correlations from lOZ equation (4) 
requires the knowledge of ^(93). The latter is obtained 
from the first member of the BGY hierarchy [26], 
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Eqs. (4), (5), and (7) constitute a very compli- 
cated set of coupled lOZ/SMSA/BGY nonhnear integro- 
differential equations with respect to h (or 5), c, and ^. 
The main problem in solving it is that the unknowns h 
and c depend on up to three variables. This leads to 
unresolvable numerical difficulties, and thus a method is 
needed to remedy such a situation. 

Any periodic function of two angle variables can be 
expanded in sine and cosine harmonics as 



f{r,ipi,ip2)= fnmH'{r)Tnl{ipi)Tml'{ip2) (8) 

n,m=0 i,;'=0,l 

using the orthogonal Chebyshev polynomials Tnoif) = 
cos{rup) and r„i(iyj) = dT„o(iy?)/diy9 = sin(niy9). Ex- 
pansion (8) can readily be applied to our two-body func- 
tions {/i, 5,c} = / with the simplification fnmU'=fnmAv 
because they are invariant with respect to the transfor- 
mation {ip\,ip2) (— y'l, — <^2) in view of the symmetry 
of Hamiltonian (1). Then exploiting the orthonormal- 
ity condition f^'^ Tni{ip)Tmi> {ip)dip = tn^nmh', where 
f„ = 7r(l — (5„o) + 27r^„0) yields the expansion coefficients 

fnml{r)=—^llf{r,(pi,<P2)Tnl{ipi)Tml{<P2)d(pidip2. (9) 

In terms of these coefficients the lOZ equation (4) re- 
duces to 



l{k)^n'm'lhn'ml{k) , (10) 



where ^nmi = ^ Jo'' i,{f)Tni{'p)Tmi{'p)<^<P are the mo- 
ments of ^(v), and the 3D Fourier transform f{k) = 
Iv /('') exp(ik • r)dr has been used. The algebraic repre- 
sentation (10) looks like the OZ equation corresponding 
to a mixture of ordinary homogeneous fluids. This is a 
very important fc^atiire because the problem can now be 
solved by adapting algorithms already known for homo- 
geneous systems. 

Furthermore, we perform the one-body polynomial ex- 
pansion 



ln^{(p) = 13 H cos Lp + Yo'nTna{f) : 



(11) 
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where only cosine harmonics appear due to the property 
^(— (p) = S,{ip)- Then the cumbersome integro-differential 
equation (7) allows to be solved analytically. 
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for n > 1, while the coefficient ao is determined from the 
normalization ^ (^{if)dip = 1. 

Handling the SMSA closure (5) also presents no diffi- 
culties, because for distances r > l^l^a (where (^(r) = 0) 
we obtain from Eqs. (5) and (6) that c(r, (^2) = 
/3[/(r)-|- J(r) cos(v3i— <^2)]- Taking into account the equal- 
ity cos(93i - 932) = Tio(¥'i)Tio(v52)+Tii((pi)Tii((p2), one 
finds cooo(r) = /9/(r-) and ciio(r) = ciii(r) = /3./(r), 
while all other c-coefficients will be equal to zero at 
r > 2^/^cr. For r < 2^/^cr, we should perform numerical 
integration (see Eq. (9)) of the right-hand side of Eq. (5) 
in order to obtain the expansion coefficients gnmi{i")- 

Another important feature is that only a small number 
Af of harmonics should be, in fact, involved because the 
expansion coefficients rapidly tend to zero with increas- 
ing jV. Then the sums X^^^ can be replaced without loss 
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of precision by finite ones with n,m < Af. In our case 
the anisotropic potential is presented by zeroth and first 
harmonics (see above the expansion for cos((pi — 1P2)), 
while a slight anharmonicity {Af > 1) in the correlation 
functions appears due to the nonlinearity of the closure. 

Once the expansion coefficients are found, all the mag- 
netic and thermodynamic properties of the system are 
obtained in a straightforward way. In particular, the 
pressure P can be calculated from the virial equation 



to critical points more closely (sec Fig. 1) and should be 
considered as more preferable than the GEMC method. 
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whereas the magnetization M = ^ J^'^ cos(ip)S_(ip)dip = 
^100- Then the phase coexistence densities between gas 
and liquid states can be evaluated by applying the well- 
known Maxwell construction to Eq. (13), while the P-F 
transition will correspond to a boundary Curie curve in 
the temperature-density plane where nonzero (sponta- 
neous) magnetization M becomes possible at H = 0. 

The coupled set of homogeneous OZ/SMSA/BGY 
equations (5), (10), and (12) was solved by adapting 
the algorithm used in Ref. [5]. The integration with re- 
spect to angle variables has been performed by Gauss- 
Chebyshev quadratures. The number of harmonics in- 
volved was Af = 3. Further increase of Af does not affect 
the solutions. Other computational details are similar 
to those of Ref. [5] when solving the Ising IE equations. 
The dimensionless quantities p* ~ pa^, T* = fceT/e, and 
H* = H/e were chosen in the presentation of the results. 

The simulations were carried out using the Gibbs 
ensemble MC (GEMC) [28] and multiple histogram 
reweighting (MHR) [29] techniques for evaluating the G- 
L and liquid-liquid (L-L) coexistences, while the Binder 
crossing scheme [24,30] was utilized to determine the P-F 
magnetic transition (at H = 0). Other simulation details 
are similar to those reported in Refs. [4,5,24]. 

In Fig. 1 we compare the OZ/SMSA/BGY results for 
the ideal (R = 00) XY fluid at different external fields H 
with the GEMC and MHR simulation data. At = 0, 
a tricritical (TC) point separates the second order P-F 
magnetic phase transition line from the first order transi- 
tion between a P-gas and an F-liquid. The iJ-dependence 
of the G-L critical temperature and density is nonmono- 
tonic. Samples of the MF binodals are included in Fig. 1 
as well to demonstrate the obvious advantage of the IE 
theory over the MF approach. Note that contrary to the 
theoretical binodals, the GEMC and MHR coexistence 
curves terminate when approaching the critical regions. 
This is because of the appearance of huge density fluctua- 
tions which cannot be properly handled within finite sim- 
ulation boxes. The MHR technique allows to approach 
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FIG. 1. The G-L binodals obtained for the ideal XY fluid 
within the OZ/SMSA/BGY approach (full curves) versus the 
GEMC (open circles) and MHR (full circles) simulation data. 

The P-F transition is plotted by the short-dashed line. The 
MF samples are shown by low lying long-dashed curves. 

The OZ/SMSA/BGY and MHR phase diagrams of the 
nonideal XY-fluid are shown in Figs. 2 and 3 for dif- 
ferent ratios R and magnetic fields H. Four types of 
phase diagram topology can be identified overall. For 
large R > 0.415 (type I), the system exhibits an ideal- 
like behavior with the existence of a TC point at H = 
and G-L transitions at H ^ for each R (Fig. 2(a)). 
At moderate values 0.26 < R < 0.415 (type II), the 
transition between a P-liquid and an F-liquid arises at 
H = additionally to the transition between a P-gas 
and a P-liquid. Here a triple point (TP) occurs too, 
where a rare P-gas, a moderately dense P-liquid, and 
a highly dense F-liquid all coexist at the same T and P 
(see Figs. 2(b) and 3(a)). The TPs can exist at H ^ 
as well and describe then the phase coexistence between 
a weakly magnetized g moderate magnetized liquid 
and a strongly magnetized hquid (Fig. 3(b)). With in- 
creasing H, either the G-L (0.376 < i? < 0.415, type Ila) 
or L-L (0.26 < R < 0.376, type lib) transition disap- 
pears in a critical end (CE) point at some finite H. For 
instance, even if R is slightly smaller than the boundary 
value i?vL = 0.376, namely R = 0.37, the L-L critical 
point ends at some H* ^ 1, while the G-L critical point 
extends to infinite field (Fig. 3(a)). In the special case 
R = Ryh, the G-L and L-L transition lines merge into 
the TC van Laar point at H* = 1.9 (Fig. 3(b)). For 
small R < 0.26 (type HI), the translational interaction 
dominates over the spin one, remaining the G-L transi- 
tion, whereas the TC point at H = transforms into a 
CE point (Fig. 2(b)). For H ^ 00, the system at any R 
behaves like a simple fluid with u{r) = (p{r) — I{r) — J{r) 
(then all spins align along H). 

As can be seen, the agreement between the theory pro- 
posed and the simulations is quite satisfactory. Slight 
deviations appear only in the vicinity of critical points. 
This is explained by finite size effects in the simulations 
and an approximate character of the SMSA closure used 
in the theory. For the latter reason, the classical value 
(3 = 1/2 of the critical exponent describing the G-L bin- 
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odal behavior |p — pd ^ \T — Tc\^ near the criticial point 
(pc) Tc) is recovered (in particular, at i? = oo and H 0), 
instead of the value /3 w 1/3 known from the renormaliza- 
tion group analysis [31]. On the other hand, the crossover 
to the TC value /3 = 1/4 can be observed near the van 
Laar point at i? = 0.376 and H* = 1.9 (Fig. 3(b)). 
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FIG. 2. The G-L and L-L binodals of the nonideal XY 
fluid within the OZ/SMSA/BGY approach (full curves) ver- 
sus the MHR data (circles). The P-F transition is plotted by 
the short- (theory) and long- (simulation) dashed curves. The 
triple point is represented by the horizontal dashed line. 
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FIG. 3. The binodals near (a) and at (b) the boundary 
value R = J?vl ■ The G-L and L-L critical points are shown in 
subset (b) for different H* as open and full squares, respec- 
tively, connected by dashed curves. The curves meet in the 
TC point (star). Other notations are the same as for Fig. 2. 

More precise IE calculations near critical points are 
possible provided a more accurate closure is used. For 
instance, the self-consistent OZ ansatz (SCOZA) [7-9] 
(which in its present formulation was implemented only 
for simple homogeneous hard-sphere Yukawa systems) 
can be extended to our inhomogcncous soft-core XY fluid 
by introducing a state dependent function K{p, T, H) 
into the SMSA closure. Then K is determined by the 
requirement of thermodynamic consistency between the 
energy and compressibility routes. In view of the inho- 
mogeneity and softness, this leads to a significant sophis- 
tication of the calculations. They go beyond the scope of 
the present Letter and will be considered elsewhere. 

In conclusion, we point out that a novel technique to 
study orientationally ordered fluids with planar spins has 
been proposed. It combines the standard IE method with 
appropriate expansions of the inhomogcncous correlation 
functions in terms of orthogonal polynomials. This re- 
duces the calculations to those inherent in a mixture of 
ordinary homogeneous fluids and thus presents now no 



numerical difficulties. Detailed comparisons with simu- 
lations have shown that the proposed approach is pow- 
erful enough to give a quantitative description of phase 
transitions in the XY spin fluid systems. 
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